/*
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#include <src/math/Operators/delta_particle.h>
#include <src/algorithm/Evolution_Algorithm/Crank_Nicolson.h>
#include <src/math/CauchyProblems/CauchyProblem_traits.h>



#if !defined(_CauchyProblem_HeatEquation_Particle_)
#define _CauchyProblem_HeatEquation_Particle_

#define NUMBER_OF_PARTICLES 20
#define DIM_D 2

template <class data, std::size_t Dim = data::Dim>
struct CrankNicolson_heat_p {
	OS_STATIC_CHECK(false); // only dim == 1 or dim == 2.	
};

template <class data>
struct CrankNicolson_heat_p<data,2> : public Crank_Nicolson<Delta_p_t<data>>
{
	typedef Delta_p_t<data>												time_operator_type;
	typedef Crank_Nicolson<Delta_p_t<data>>									Base;

	virtual OS_double CFL() {
		return controled_moment_cfl<time_operator_type>()(Op_iterate);
	};
	virtual bool criteria() {
		OS_double breakpoint1 = jmmath::norm22(Op_iterate - Op_S_k1);
		OS_double breakpoint2 = jmmath::norm22(Op_iterate);
		return breakpoint1/(1.+breakpoint2) < 1e-4;
	};
};

template <class data>
struct CrankNicolson_heat_p<data,1> : public Crank_Nicolson<Delta_p_t<data>>
{
	typedef Delta_p_t<data>												time_operator_type;
	typedef Crank_Nicolson<Delta_p_t<data>>									Base;

	virtual OS_double CFL() {
		return controled_moment_cfl<time_operator_type>()(Op_iterate);
	};
	virtual bool criteria() {
		OS_double breakpoint1 = jmmath::norm22(Op_iterate - Op_S_k1);
		OS_double breakpoint2 = jmmath::norm22(Op_iterate);
		return breakpoint1/(1.+breakpoint2) < 1e-4;
	};
};


template <class data, class Interface = Graphical_Interface>
class CauchyProblem_HeatEquation_traits : public CauchyProblem_traits<  CrankNicolson_heat_p <data>, Interface>
{
public :
	typedef CauchyProblem_traits< Crank_Nicolson < Delta_p_t<data> > , Interface> Base;
	CauchyProblem_HeatEquation_traits() : dim_d_(DIM_D)
	{
	};

protected :

	OS_int dim_d_;
};

/*!	\brief Multidimensionnel Cauchy Problem \f$\partial_t u(t,x) = \Delta u(t,x)\f$.
	This class initialize the Cauchy problem for the multi dimensional heat equation
	\f[
		\partial_t u(t,x) = \Delta u(t,x)
	\f];
	with \f$u(0,x) = u_0(x) > 0\f$, \f$\int_{\RR^d}u_0 = 1\f$, \f$x \in \RR^d\f$, \f$t\ge 0\f$.

	The numerical method is a particle one (Optimal scheme) : it consists in approximating the optimal quantizers
	by the following dynamical system :
	\f[
		 \frac{d}{dt} S_i = d \sum_{j \in \mathcal{V}_{i} }  \left(\omega_{i,j}\right) \frac{S_i-S_j}{\|S_i-S_j\|^2}
		 S_i(0) = S_i^0.
	\f];

 with the following notations and assumptions :
 -
 - \f$\{\frac{1}{N},S_i(0)\}_{i=1,..,N}\f$ is an optimal quantizer of the initial density probability \f$u_0\f$. Thus we assume
 that the initial condition admits an equi probable quantizer.
  - \f$\mathcal{V}\f$ is the Delaunay graph induced by \f$\{S_i\}_{i=1,..,N}\f$. \f$j \in \mathcal{V}_{i}\f$ are all indices of particles connected to \f$S_i\f$ in the Delaunay graph.
\f$\omega_{i,j}\f$ is the weight of each vertex  : we choose to define
\f$\omega_{i,j} = \frac{1}{|\mathcal{V}_{i}|})\f$, where \f$|\mathcal{V}_{i}|\f$, the cardinal of \f$\mathcal{V}_{i}\f$, holds
for the number of particles connected to \f$S_i\f$. 
- \f$\|S_i-S_j\|\f$ holds for the standard euclidian norm.

This dynamical system is numerically approximated by a standard Crank Nicolson scheme :
	\f[
		\delta_t S_i^n = d \sum_{j \in \mathcal{V}_{i}^n }  \left(\omega_{i,j}^n\right) \frac{S_i^{n+1/2}-S_j^{n+1/2}}{\|S_i^{n+1/2}-S_j^{n+1/2}\|^2}
	\f];

 
where \f$\delta_t S_i^n := \frac{S_i^{n+1}-S_i^n}{t^{n+1}-t^n} and S_i^{n+1/2} = \frac{S_i^{n+1}+S_i^{n}}{2}\f$

 The initial data \f$\{S_i^n\}_{i=1,..,N}\f$ corresponds to a normally distributed initial conditions. More precisely,
 it is computed by the Box Muller method together with a Sobol generator.
*/
		template <class data, std::size_t Dim = data::Dim, class Interface = Graphical_Interface>
		class Cauchy_Delta_p : public CauchyProblem_HeatEquation_traits< data, Graphical_Interface>
		{
			OS_STATIC_CHECK(data::Dim > 1);

			public :
			virtual smart_ptr<interface_type> Create_Instance() {
				return smart_ptr<interface_type>(new Cauchy_Delta_p<data,Dim>);
			};

			Cauchy_Delta_p(){
				initial_conditions(NUMBER_OF_PARTICLES);
			};
			virtual void set_dim(OS_size dim_d)
			{
				OS_DYNAMIC_CHECK(dim_d > 1,"Only for dimension greater or equal than 2");
				dim_d_ = dim_d;
				initial_conditions(NUMBER_OF_PARTICLES);
			};
			OS_size get_dim() {return dim_d_;};
			void initial_conditions(OS_size nb_part)
			{
				set_initial_conditions(InverseCumulativeFunctor<data>()(nb_part,dim_d_));
			};
		};

/*! \brief One dimensional Cauchy Problem - Particle method
    This class initialize the Cauchy problem for the one dimensional heat equation
\f[ \label{heat_e}
    \partial_t u = \partial_x^2 u; u(0,x) = u_0(x) > 0, \int_{\RR} u_0 = 1
\f]
    with \f$x \in \RR, t\ge 0\f$.

 We provide a one dimensional particle method of this heat equation.
 We recall that the particle method for the heat equation consists in considering the equivalent equation
\f[ \label{heat_p}
    \partial_t S = - \partial_y {\frac{1}{\partial_y S}}; S(0,y) = S_0(y)
\f]
with \f$y \in [0,1]\f$. \f$S(t,y)\f$ is here the inverse of the cumulative of \f$u(t,x)\f$, that we denote \f$S^{-1}(t,x)=\int_{-infty}^xu(t,s)ds\f$.
The equation (\ref{heat_p}) is approximated with the following
 dynamical system of \f$N\f$ particles \f$ \{S_i \}_{i=1,..,N} \f$ :
\f[ \label{sd_heat}
    \frac{d}{dt} S_i = - \delta_{\bar{y}} {\frac{1}{\delta_{y} S_i}}; S_i(0) = S_i^0
\f]
 where \f$\{S_i \equiv S_i(t)\}_{i=1,..,N}\f$ are the particles location,
\f$ \delta_{y} S_i = \frac{S_{i+1}-S_i}{h}\f$, \f$\delta_{\bar{y}} S_i = \frac{S_{i}-S_{i-1}}{h}\f$
are the (normalized) difference operators, with \f$h:= \frac{1}{N}\f$.

  \f$\{S_i^0\}_{i=1,..,N}\f$ is the optimal quantizer of \f$u_0\f$ (we assume that
 the optimal quantizer uf \f$u_0\f$ is equiprobable) : i.e. \f$S_i^0 = S^0 { \frac{i+1/2}{N}}, i=0,..,N-1\f$,
where \f$S^0\f$ is the inverse function of the cumulative of \f$u_0\f$, that we denote \f$(S^0)^{-1}(x) = \int_{-\infty}^x u_0(s)ds\f$.

 The hamilton equations (\ref{sd_heat}) is approximated directly by the non linear system of equations
 \f[
    \delta_t S_i^n := \frac{S_i^{n+1}-S_i^n}{\tau^{n}} = - \delta_{\bar{y}} {\frac{1}{\delta_{y} S_i^{n+1/2}}}
\f]
 where \f$S_i^{n} \sim S_i(t^{n})\f$ and  \f$S_i^{n+1/2} = \frac{S_i^n+S_i^{n+1}}{2}\f$. We know that for \f$\tau^n = t^{n+1}-t^n\f$ small enough, the previous equation has
 an unique solution. The choice of Crank Nicolson scheme is motivated by the fact that it will verify
the entropy evolution property
 \f[
    \delta_t \frac{1}{2N}\sum_i |S_i^n|^2 = 1
\f]
that corresponds to the property of the original equations
 \f[
    \frac{d}{dt} \int_{[0,1]} S^2(y) dy = 1 = \int_{\RR} x^2 u(t,x) dx
\f]
Thus this scheme suits well to our purposes.


Here, the initial condition correspond to a cauchy data given by a concentrated gaussian
\f$u_0(x) = N(x)_{0,t} = \frac{1}{\sqrt{2 t}} \exp {\frac{-x^2}{t}}\f$, where t=0.1.
This leads to initial condition given by
\f$\{S_i^0 = S^{-1}_{0}(i+1/2)\}_i\f$, where \f$S^{-1}_0(x)\f$ is the cumulative of the gaussian.
*/

		template <class data, class Interface>
		class Cauchy_Delta_p<data,1,Interface> : public CauchyProblem_HeatEquation_traits<data,Interface>
		{
		public :
			OS_STATIC_CHECK(data::Dim == 1);

			typedef Cauchy_Delta_p< data,1> cauchy_problem_type;
			virtual smart_ptr<interface_type> Create_Instance() {
				return smart_ptr<interface_type>(new cauchy_problem_type);
			};
			virtual OS_size get_dim() {return 1;};


			Cauchy_Delta_p(){
				initial_conditions(NUMBER_OF_PARTICLES);
			};
			void set_dim(OS_size dim) {
				OS_DYNAMIC_CHECK(dim == 1,"Cauchy_Delta_p<data,1> : only dimension 1");
			};

			virtual void initial_conditions(OS_size nb_part)
			{
				set_initial_conditions(InverseCumulativeFunctor<data>()(nb_part));
			};
		};
#undef NUMBER_OF_PARTICLES
#undef DIM_D
#endif